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Black holes binaries support unstable orbits at very close separations. In the simplest case of 
geodesies around a Schwarzschild black hole the orbits, though unstable, are regular. Under per- 
turbation the unstable orbits can become the locus of chaos. All unstable orbits, whether regular 
or chaotic, can be quantified by their Lyapunov exponents. The exponents are observationally rel- 
evant since the phase of gravitational waves can decohere in a Lyapunov time. If the timescale for 
dissipation due to gravitational waves is shorter than the Lyapunov time, chaos will be damped and 
essentially unobservable. We find the two timescales can be comparable. We emphasize that the 
Lyapunov exponents must only be used cautiously for several reasons: they are relative and depend 
on the coordinate system used, they vary from orbit to orbit, and finally they can be deceptively 
diluted by transient behaviour for orbits which pass in and out of unstable regions. 
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Newtonian gravity predicts the elliptical planetary or- 
bits around the sun which Kepler described. Einstein 
gravity predicts precessing elliptical orbits around a cen- 
tral star, thereby reconciling the precession of the peri- 
helion of Mercury with general relativity. In the extreme 
case of a central black hole, there are also a simple set 
of unstable circular orbits in addition to the usual sta- 
ble circular orbits. Related to these are the homoclinic 
orbits which lie on the boundary between dynamical sta- 
bility and instability [1,2]. To the list of possible orbits, a 
set of chaotic orbits has recently been added for rapidly 
spinning black holes [3-5] (and for the interesting but 
less physically realistic Majumdar-Papapetrou black hole 
pairs of equal mass and charge [6-9]). With the future 
gravitational wave experiments LIGO and LISA we hope 
to see these innermost orbits and reconstruct a map of 
spacetime around black holes. 

The instability of the innermost orbits around black 
holes will etch certain landmarks in a gravitational wave 
map. Gravitational waveforms of neighboring orbits will 
dechohere in a time scale set by the instability [10]. If the 
timescale for dissipation through gravitational radiation 
is faster than the instability timescale, then chaos will 
be damped and the gravitational wave signal will not 
observably decohere. 

The simple set of unstable circular orbits around a 
Schwarzschild black hole are a consequence of the non- 
linearity of general relativity. Their instability can be 
quantified by a positive Lyapunov exponent [10]. Al- 
though Lyapunov exponents are often associated with 
chaotic dynamics, the geodesies around a Schwarzschild 
black hole are not chaotic: the orbits are fully soluble 
and therefore integrable. However, under perturbation, 
chaos is likely to develop along the unstable circular and 
homoclinic orbits. An example of this has been found 
when the black holes spin. The nonlinearity degener- 
ates to a nonintegrability and chaos [3-5]. The number 
of unstable periodic orbits proliferates so that they have 
to pack themselves into a fractal in order to crowd into 



that region of phase [7,5]. The unstable orbits will have 
positive Lyapunov exponents [11,12] and will emerge as 
fractals in phase space [4,5]. 

The Lyapunov exponents, while a seemingly useful 
tool, have uncomfortable shortcomings in the context of 
general relativity. Firstly, the Lyapunov exponents vary 
from orbit to orbit and so do not have the surveying 
power to scan the collective behaviour of all orbits that 
fractals methods do. Secondly, the Lyapunov exponents 
are a measure of the deviation of two neighboring orbits 
in time and therefore overtly depend on the time coordi- 
nate used. Since time is relative so too are the Lyapunov 
exponents. The relativity of the Lyapunov exponents has 
been known to erroneously lead to zero Lyapunov ex- 
ponents for truly chaotic systems [13-16]. Importantly, 
topological measures of chaos such as fractals are coordi- 
nate invariant and are not plagued by the relativism of 
space and time [7,15]. 

In rare cases when there is a prefered time direction the 
ambiguity of time can be avoided. For the simplest case 
of a Schwarzschild black hole there is a timelike Killing 
vector which selects a prefered time direction. In other 
words, from our position asymptotically far away from 
the black hole, we use a well defined time coordinate in 
our observations. As long as we conscientiously compare 
all timescales in the same coordinate system, we should 
get meaningful comparisons. 

We investigate the stability of three types of black hole 
binary: (i) Schwarzschild black hole (non-spinning, test 
particle motion), (ii) the Post-Newtonian (PN) expansion 
of the two-body problem (non-spinning black holes), and 
(iii) the chaotic orbits of spinning black holes in the PN- 
expansion. 

The stability analysis begins with the equations of mo- 
tion summarized as 



(1.1) 



To analyze the stability of a given orbit we linearize the 
equations of motion about that orbit 
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with 
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Xi(t) 



the linear stability matrix. The solution to the linearized 
equations can be written as 

5X i {t)=L ij {t)8X j {Q) (1.4) 

in terms of the evolution matrix which must obey 

Lij(t) — Ki rn L m jit) (1-5) 



and Lij(0) = dij. A determination of the eigenvalues of 
Li 

cally 



j t] leads to the principal Lyapunov exponent. Specifi- 
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II. SCHWARZSCHILD ORBITS 
A. Circular orbits 

Here we evaluate the Lyapunov exponent of unstable 
orbits around the Schwarzschild black hole. We consider 
the usual geodesies of a non-spinning, light companion. 
We work in Schwarzschild time, the time measured by 
an observer asymptotically far from the black hole. The 
Lyapunov exponent was already evaluated in Ref. [10] 
in a different time coordinate system. This exemplifies 
the ambiguity of time. Still, the timescales which were 
compared in that paper were all measured in the same 
coordinate system and therefore the general conclusions 
of Ref. [10] still hold. 

To isolate A, we begin with the Lagrangian for a (non- 
spinning) test particle in the Schwarzschild spacetime 
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The black hole mass has been set to unity. We consider 
motion in an equatorial plane to eliminate the cylic 9 
variable by setting it equal to it/ 2 and define the canon- 
ical momenta by 8C/5idq/ds) = p q : 
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To change into Schwarzschild time t, we use eqn. (2.2) to 
define the transformation 
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Notice that in Ref. [10] an unusual time coordinate t' was 
used instead which was defined by the transformation 
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We will redo the stability analysis in Schwarzschild time 
t. The equations of motion can be derived through 



Sdr/ds Sr 
and reduce to a two-dimensional system: 
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where an overdot denotes differentiation with respect to 
Schwarzschild time t. To compare with Ref. [10], we first 
consider circular orbits. Linearizing the equations of mo- 
tion with Xi{t) — (p r , r) about orbits of constant r gives 
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The eigenvalues along circular orbits are 
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For the unstable circular orbit at r = 4 the angular 
momentum and energy are L = 4, E = 1 respectively 
and the eigenvalues of Kij are 



A± =± 
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(2.9) 



The conservation of energy ensures that in these canon- 
ical coordinates, the Lyapunov exponents must come in 
plus-minus pairs to conserve the volume of phase space. 

The unit normalized eigenvectors corresponding to A ± 
are 



e + = i(-l,2V2) 
e_ = i(l,2V2). 



In this eigenbasis K t j is diagonal with 
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so that 
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Notice that in Ref. [10], the analysis was carried out in a 



four-dimensional coordinate system: X[ = (p r 



,0)- 



The additional coordinates are unimportant in the dy- 
namical study and only hampered the diagonalization of 
Kij [10]. If we were to redo the stability analysis in 
three coordinates X" = (p r ,P<f,, r), we would get the same 
eigenvalues (2.9) and an additional A = 0. If we move up 
to the four coordinates of X ■ we add yet another A = 
giving a degenerate set of eigenvalues. The matrix 
cannot be diagonalized in the event of degenerate eigen- 
values which leads to an unnecessary complication. For 
this reason, we stick to the pertinent two-dimensional 
system Xi = (p r , r). 

The relativity of time and the influence on the 
Lyapunov exponent is apparent at this stage. In 
Schwarschild time t we find A± = ±^75 while in the 
time t' used in Ref. [10], at r — 4 the exponents were 
found to be A' ± = The Lyapunov exponents are 



not coordinate invariant. However at r 
it follows that the combination 



At = X't' 



4, t' = t/A and 
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is invariant. 

We compare the Lyapunov timescale T\ = 1 / A to the 
gravitational wave timescale T w = 2ir/(j). For the or- 
bit at r = 4, T\/T w = y/2/n w 0.45. The Lyapunov 
timescale is less than about one orbit around the cen- 
tral black hole. Notice that even though Ref. [10] oper- 
ated in an unusual time, the calculations were performed 
self-consistently so that the ratio of T\ to T w is correct. 
The Lyapunov timescale is shorter than the gravitational 
wave timescale, making the instability observationally 
relevant. 

In principle we could also compare T\ to the de- 
cay time due to energy lost in the form of gravita- 
tional radiation. For a test-particle in a circular or- 
bit around a Schwarzschild black hole the decay time is 
Td = (5/256)r 4 //t where /j, is the reduced mass. At r = 4 
this is = 5/ (3/i ) 3> 1 since in the test-particle limit 
and T\ will be shorter than the decay time, again 
making the instability observationally relevant. 



B. Homoclinic orbits 



systems point of view, the true division between stability 
and instability occurs more generally along the homo- 
clinic orbits [1,2]. 

The stable manifold of a periodic orbit is defined as 
the set of points in phase space that, when evolved for- 
wards in time, approach the periodic orbit. The unstable 
manifold is the set of points in phase that, when evolved 
backwards in time, approach the periodic orbit. For an 
integrable, nonchaotic system, the stable and unstable 
manifolds can intersect each other along a single orbit. 
This orbit is called homoclinic if it approaches the same 
fixed point in the past and in the future. For black holes 
the homoclinic orbits begin at an unstable circular orbit, 
roll out to a maximum radius and fall back in to the same 
unstable circular orbit. An example is shown in fig. 1. 
The homoclinic orbits are sometimes called zoom-and- 
whirl orbits in the gravitational wave literature because 
the orbits whirl around the center of mass and then zoom 
out into an ellipse before whirling in again. 




FIG. 1. A segment of the homoclinic orbit with j3 = 1/4. 

Under perturbation the homoclinic orbits can become 
the site of a homoclinic tangle. The tangle occurs when 
the stable and unstable manifolds intersect transversely 
at an infinite number of points. The intersection will no 
longer occur along a simple line but will instead define a 
fractal set of chaotic orbits. In this section, we study the 
stability of the nonchaotic, simple set of homoclinic or- 
bits around a Schwarzschild black hole. In §IV we study 
chaotic orbits of spinning black holes. 

As is well known, orbital motion around a 
Schwarzschild black hole conveniently reduces to one- 
dimensional motion in an effective potential 



Although the unstable circular obits are often em- 
phasized they are actually a subset of the pertinent or- 
bits. The division between stability and instability for a 
Schwarzschild black hole is often taken to be the inner- 
most stable circular orbit (ISCO). The ISCO is actually 
the saddle point at which the unstable circular orbit co- 
incides with the stable circular orbit. From a dynamical 
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Circular orbits are solutions of V e s = E. For large 
enough angular momentum there are two circular orbits, 
one unstable and one stable. As found in Ref. [1], the 
homoclinic orbits have E < 1 and are described by the 
solution 
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where < (3 < 1/2 and t{<p) is a complicated function 
[1]. * The circular orbits can also be parameterzied by (3 



as 



^unstable = 6/(1 + /3) 
^stable =6/(1-/3) 



L = 2y3/ (1 — /3 2 
2-/3 



v/2/(l-/3). 




FIG. 2. The effective potential for orbits with /3 = 1/4. 
The segment in bold marks the corresponding homoclinic tra- 
jectory. The homoclinic orbit begins at the unstable circular 
orbit at the top of the hill (r = 4.8), rolls out past the stable 
circular orbit in the valley and on out to the maximum radius 
(r = 12) before rolling in again and climbing back up to the 
unstable orbit. 

A homoclinic oribt starts at r uns t a bio and rolls out to 

W - 6/(1 - 2/3) (2.17) 

winding around the black hole as it does so and then 
drops back in to r unsta bi c . The ISCO is a homoclinic 



"There appears to be typo in eqn. (1.5) of [2]. Eqn. (2.16) 
has a factor 1 /2 which is missing from the second term in eqn. 
(1.5). 



orbit with (3 = 0. The first homoclinic orbit at (3 = 1/2 
starts at the unstable orbit at r = 4 and rolls out to 
infinity before returning. For (3 = 1/4, r unsta bio = 4.8, 
f stable = 8 and r max = 12. The effective potential for 
the homoclinic orbit is drawn in fig. 2. A segment of the 
orbit in the equatorial plane is represented in fig. 1. 
To analyze the stability we linearize to find 
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The most general eigenvalues are 
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Strictly speaking, I is a stability exponent and is not iden- 
tical to the time averaged Lyapunov exponent defined 
in eqn. (1.6). Figure 3 shows the real and imaginary 
parts of the positive stability exponent. As expected, 
the exponent is postive near the unstable inner radius 
and becomes imaginary in the vicinity of the stable cir- 
cular radius dropping down to nearly zero as it reaches 
the apihelion and then runs back through these values as 
it moves back in to perihelion. 



FIG. 3. 13 = 1/4 The solid line is the Real part of the pos- 
itive stability exponent and the dotted line is the Imaginary 
part. 

Because of this time variability in the stability along 
the orbit, we have to be cautious in interpreting the 
timescales. The gravitational wave frequency will jag up 
and down as the orbit zooms and whirls [2]. And it isn't 
obvious which timescales to compare. Instead of using 
the variable, analytic result we could try a time average. 

To this end we compare the analytic value of I from 
eqn. (2.19) to the time average Lyapunov exponent de- 
fined from eqn. (1.6) 
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For comparison we first look along the unstable circu- 
lar orbit at r — 4.8. The upper panel of figure 4 shows 
the Lyapunov exponent as defined by eqn. (1.6) for the 
unstable circular orbit. The exponent was obtained by 
numerically integrating the Lij using eqn. (1.5) and fig- 
ure 4 shows Xt versus t from eqn. (1.6). The numerically 
calculated value shown in figure 4 is identical to the an- 
alytic value given by eqn. (2.8) of A m 0.475. 

However for a homoclinic orbit which begins at r = 4.8, 
the time averaged (1.6) behaves as though there is no 
instability (lower panel of figure 4), when we know from 
the analytic result shown in fig. 3 that there is. The 
time averaged Lyapunov exponent will vanish along this 
orbit even though it clearly has an unstable segment. We 
have to be very cautious therefore when we interpret the 
Lyapunov exponent. 
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FIG. 4. (3 = 1/4. The slope of the line in the top panel is 
the numerically calculated Lyapunov exponent A ~ 0.475 for 
the unstable circular orbit at r — 4.8 using eqn. (1.6). This 
matches exactly the analytic value predicted from eqn. (2.8). 
The slope of the line in the lower panel shows zero Lyapunov 
exponent as calculated by eqn. (1.6) for the homoclinic orbit. 
This is to be contrasted with the analytic value of the stability 
exponent shown in fig. 3 which has positive segments 

For emphasis, if one was just scanning numerically, 
the mistaken conclusion could be drawn that these or- 
bits were dynamically simple. This may turn out to be 
particularly important for the chaotic orbits of §IV. 



III. POST-NEWTONIAN ORBITS 

To move beyond the test particle limit, the two-body 
problem has been expanded in a Post-Newtonian (PN) 
expansion approximation to the fully relativistic two- 
body problem [17-19]. In this section we consider two 
black holes which are not spinning. In Ref. [2] the stabil- 
ity of the fixed points in the PN equations to second-order 
(2PN) was tested following [17]. We quote the results of 
Ref. [2] here. To second-order in the PN expansion, the 
center of mass equations of motion for the binary orbit 
can be written in harmonic coordinates as [17-19] 



rh = r h (, 



±(A + Br h ) 



\B + 2^ 



(3-1) 
(3.2) 



where M = 1 is the total mass of the pair. The transfor- 
mation between harmonic coordinates and Schwarzschild 
coordinates is = r — m. The form of A(rh, fh, <f) and 
B(rh,fh,<j>) depends on the relative masses of the two 
black holes and on the order of the PN expansion and 
can be found in Ref. [17]. 

As in Ref. [17], the stability of the fixed points is tested 
by perturbing eqns. (3.1)-(3.2) about a circular orbit to 
obtain, 
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where r^ Q is the radius of the circular orbit in harmonic 
coordinates and <j)g = mA jr\ is a function of the ra- 
dius of the orbit and is found explicitly in Ref. [2]. The 
eigenvalues of (3.3) are, 



t = 0, 



±(a+ be) 



1/2 



(3.5) 



Stable oscillations about a circular orbit correspond to 
imaginary I so that a + be < 0. Unstable orbits corre- 
spond to real positive t and so have a + bc > 0. The value 
of £ 2 for equal mass binaries as a function of the circular 
radius is plotted in fig. 5. 
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FIG. 5. £ 2 as a function of the constant circular radius 
in harmonic coordinates, r^o- i 2 > corresponds to unstable 
circular orbits and I 2 < corresponds to stable circular orbits. 



Using the results of Ref. [2] we deduce that in the test- 
mass limit, the innermost unstable circular orbit occurs at 
Schwarzschild radius r a = r] lo + 1 = 4.96. A comparison 
to the gravitational wave timescale in the PN expansion 
then gives T\/T w w 0.158. In the opposite extreme of 
equal mass binaries, the innermost unstable circular orbit 
occurs at r Q = r ho + l = 5.78 and T x /T w w 0.21. A direct 
comparison to the Schwarzschild case isn't that valuable. 
What is noteworthy is that the Lyapunov timescales are 
again less than about one orbit around the center of mass. 
Consequently, one expects the decoherence of the gravi- 
tational waveform to be observationally significant. The 
instability timescales are comparable to the decay times 
although of course even a small loss in energy can induce 
merger for such an unstable orbit. 

The analysis extended to homoclinic orbits can be 
gleaned from Ref. [2]. The homoclinic orbits to 2PN or- 
der show similar features to the homoclinic orbits of the 
Schwarzschild spacetime. The analytic Lyapunov expo- 
nent will pass from positive to imaginary values as the 
orbit winds around the center of mass. However, the 
time-averaged exponent will dilute these critical features. 

None of these orbits are chaotic although they are un- 
stable. We turn to the chaotic orbits of rapidly spinning 
binaries next. 



IV. CHAOTIC ORBITS 

The dynamics can become chaotic when the homoclinic 
orbit is perturbed leading to a homoclinic tangle. The 
intersection of the stable and unstable manifold will no 
longer occur along a line but will intersect transversally 
an infinite number of times. The fractal set of unstable 
chaotic orbits lies along this tangled intersection. Ref. [1] 
studied generic gravitational perturbations along the ho- 
moclinic orbits and found, as they expected, that the dy- 
namics could become chaotic. The physical significance 
of the perturbations however wasn't clear and therefore 
the observational consequences were difficult to assess. 




2x10 4 4x10" 
FIG. 6. The slope of the line is the Lyapunov exponent for 
the chaotic orbit of fig. 7. 

The chaotic dynamics discovered in Ref. [3] for a supra- 
maximally spinning test-particle and for rapidly spinning 
pairs in the Post-Newtonian expansion [4,5], may occur 
along these homoclinic orbits. At the least the chaotic 
behaviour kicks up most conspicuously in the vicinity of 
these orbits. 

We cannot determine the Lyapunov exponents analyt- 
ically since the orbits are not analytically soluble - the 
very meaning of nonintegrability. We can however use 
eqn. (1.6) to numerically determine the exponents as we 
have done in Ref. [12]. The Lyapunov exponent for a 
maximally spinning pair of black holes is shown in fig. 
6. The value read from this is T\ w 11 in units of wind- 
ings around the center of mass. A segment of the orbit 
projected onto a plane is shown in fig. 7. 




FIG. 7. A projection onto the plane of two maximally spin- 
ning black holes in a chaotic orbit. 

For such an erratic orbit, it is not simple to define 
the gravitational wave timescale or the radiation reac- 
tion timescale. By starting the numerical simulation at 
a radius greater than 20 we found a rough estimate of 
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Td ~ 4 — 5 orbits so that T\/Td > 2 and the Lyapunov 
timescale is longer than the dissipation timescale. For 
this orbit the gravitational waveform will not have time 
to decohere before plunge and observations will not be 
severly disrupted by irregularity. However a factor of 2 is 
a tight margin especially since the PN-expansion is being 
pushed to extremes at such small separations. 

It is important to stress that this is but one specific 
orbit and the ratio of T\/T w will vary from orbit to or- 
bit. Some chaotic orbits will undoubtedly have a longer 
Lyapunov timescale to dissipation timescale. A broad 
survey of the irregular region of phase space would be 
valuable. At 2PN order this is not so useful since the 
PN-expansion converges very slowly to the full relativis- 
tic problem. Notice for instance that T\ is longer for this 
two-body approximation than it is for the fully relativis- 
tic Schwarzschild orbits. This may be a consequence of 
the slow convergence of the expansion and may show an 
underestimate of the instability. It is also possible that, 
like the homoclinic orbits, the measure of instability is 
diluted by the time average across the orbit. A survey at 
higher orders may be useful but requires a greater than 
3PN-expansion that includes all spin arrangements and 
is not restricted to circular or quasi-circular orbits. We 
hope these higher orders will be available imminently and 
a survey of orbits will be viable in the near future. 

SUMMARY 

The Lyapunov exponents can be very useful for com- 
parisons of physical scales. However, Lyapunov expo- 
nents also have some shortcomings which require we 
tread cautiously: 

• Lyapunov exponents are relative. They depend on the 
worldline of the observer and the time they measure. 

• They vary from orbit to orbit and may not contain 
generic information. 

• They can give zero when averaged over orbits which 
move in and out of unstable regions. 

Therefore, while important and useful, the Lyapunov ex- 
ponents can be misleading and can only be used cau- 
tiously. 

As far as we can trust them, the Lyapunov exponents 
give an estimate of the importance of instability to obser- 
vations of gravitational waves. If the Lyapunov timescale 
is short compared to the inverse frequency of the gravita- 
tional waves emitted and is short compared to the dissi- 
pation timescale then instability will cause an observable 
decoherence of gravitational waves. We found that 

• the Lyapunov timescale is shorter than both the gravi- 
tational wave timescale and the dissipation timescale for 
unstable circular orbits in the approximation of a test- 
particle around a Schwarzschild black hole, 

• the Lyapunov timescale is shorter than the gravita- 
tional wave timesecale and comparable to the dissipation 



timescale for unstable circular orbits in the 2PN approx- 
imation in the absence of spins, 

• and that the Lyapunov time was about a factor of 2 
larger than the decay time for one randomly sampled 
chaotic orbit of a pair of maximally spinning black holes 
in the 2PN expansion. 

The longer Lyapunov time for orbits in the 2PN ap- 
proximation versus the test-particle approximation may 
be a real effect or it may be due to the slow convergence 
of the 2PN expansion to the full nonlinear problem or 
finally it may be due to the time average over such a var- 
ied orbit. In short, dissipation due to gravitational waves 
does abate chaos although the competition between chaos 
and dissipation is close. Better approximations to the 
two-body problem are needed to determine conclusively 
if chaos will affect observations of gravitational waves. 
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